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In the authors’ previous studies [1], a time-accurate, upwind finite volume method (ETAU scheme) for 
computing compressible flows on unstructured grids was proposed. The scheme is second order accurate in 
space and time and yields high resolution in the presence of discontinuities. The scheme features a multidi- 
mensional limiter and multidimensional numerical dissipation. These help to stabilize the numerical process 
and to overcome the annoying ‘pathological behaviors’ of upwind schemes. 

In the present paper, it will be further shown that such multidimensional treatments also lead to a nearly 
‘all-speed’ or ‘Mach number insensitive’ upwind scheme. For flows at very high Mach number, e.g., 10, lo- 
cal numerical instabilities or the ‘pathological behaviors’ are suppressed, while for flows at very low Mach 
number, e.g., 0.02, computation can be directly carried out without invoking preconditioning. For flows in 
different Mach number regimes, i.e., low, medium, and high Mach numbers, one only needs to adjust one or 
two parameters in the scheme. 

Several examples with low and high Mach numbers are demonstrated in this paper. Thus, the ETAU scheme 
is applicable to a broad spectrum of flow regimes ranging from high supersonic to low subsonic, appropriate 
for both CFD (computational fluid dynamics) and CAA (computational aeroacoustics). 

I. Introduction 

Upwind finite volume (FV) schemes are broadly employed in computational fluid dynamics (CFD) and computa- 
tional aeroacoustics (CAA) due to their robustness and geometric flexibility. The upwind FV schemes are based on 
the Gauss divergence theorem applied to a control volume (CV) using a Riemann solver to generate surface fluxes. 

The Godunov scheme and the TVD high resolution schemes [2, 3] are the fundamental upwind schemes. In past 
decades significant improvement of the accuracy of the upwind methods has been achieved by using higher order 
approximations, such as the ENO (essentially non-oscillatory) [4] and WENO (weighted ENO) [5] schemes, the 
DG (discontinuous Galerkin) [6] scheme, the SV (spectral finite volume) and the SD (spectral difference) schemes 
[7, 8]. Recently, the authors [1] proposed an ETAU (enhanced time accurate upwind) scheme which avoids high 
computer costs in practical computations but retains reasonable accuracy (nominally second order in space and time) 
and robustness. 

For any scheme, it is desirable to cover the entire Mach number range from extremely low Mach number to 
extremely high Mach number, without special treatment. Such so-called ‘all speed’ or ‘Mach number insensitive’ 
scheme has been a research topic in CFD for decades. For high Mach number flows, due to some inherent local 
instability, an upwind scheme may exhibit the so-called ‘pathological behaviors’, such as the carbuncle phenomenon, 
the expansion shock, etc., which spoil the numerical solutions. Numerical dissipation is required to stabilize the 
scheme. [11-14]. For extremely low Mach number flows, an upwind scheme often needs preconditioning treatment 
for a successful computation. 

When the ETAU scheme was constructed, it was designed to overcome the pathological behaviors that often occur 
at high Mach number [1]. This is achieved by applying multidimensional dissipation to the scheme and systematically 
removing such annoying phenomena [1, 15]. For different flow regimes (i.e., over a Mach number range), one only 
needs to adjust a weighing factor, p, which controls the numerical dissipation. 

Recently, further investigation of the ETAU scheme shows that the scheme is also capable of computing flows at 
very low Mach number (e.g., 0.02) without invoking preconditioning. This renders the ETAU scheme a nearly ‘all 
speed’ or ‘Mach number insensitive’ upwind scheme. Compared to the conventional upwind schemes, analysis shows 
that such capability owes to the following built-in numerical procedures in the ETAU scheme: 
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• multidimensional limiting, 

• multidimensional dissipation, 

• careful surface flux evaluation. 

With a selection of merely one or two parameters for different flow regimes, i.e., low, medium, and high Mach number, 
the ETAU scheme is capable of covering a wide range of Mach number (e.g., from 0.02 to 10), including both weak 
acoustic waves and strong shocks. 

The major features of the ETAU scheme are reviewed in Section II. Numerical examples with very low Mach 
number to very high Mach number are presented in Sections III and IV. Concluding remarks are addressed in Section 
V. 


II. The ETAU Scheme 

The enhanced time- accurate upwind (ETAU) scheme [1] is nominally second order accurate in space and time 
(in the absence of flow discontinuities). As mentioned in §1, the scheme includes many attributes often found in the 
construction of high order upwind schemes. Details about the scheme can be found in [1], but the important features 
are described below, beginning with the conservation form of the Euler/Navier-Stokes (E/N-S) equations. 

A. Conservation form of the compressible Euler/Navier-Stokes equations 

Let p, «, v, /?, and y be the density, the two velocity components, static pressure, and constant specific heat ratio, 
respectively. The two-dimensional unsteady E/N-S equations are written in the standard conservation form: 

Uf+Fx + G-y = 0, (1) 

where x, y , and t are the streamwise and transverse coordinates and time, respectively. The flux vectors are further 
split into inviscid and viscous fluxes: 


F — Fj — F v , G — Gj — Gy. 

The conservative flow variable vector U, and the inviscid flux vectors F* and G* are given in nondimensional form as: 
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Here e = + 1/2(m 2 + v 2 ), and the enthalpy is H = p/p + e. The viscous flux vectors F v and G v can be found in 
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Figure 1 . (a) A typical unstructured triangular grid in two-dimensional space, (b) control volume in £ 3 , (c) a multidimensional dissipation model. 


2 of 11 


American Institute of Aeronautics and Astronautics 


B. Cell surface flix evaluation 


A typical two-dimensional unstructured triangular cell used in the scheme is illustrated in Fig. 1. Here, A 4BC is 
a triangular cell centered at O and D,E,F are the centers of the neighboring triangular cells. The flow variables 
at the previous time step are stored at these triangle cell centers. By considering (x,y, t) as coordinates of a three- 
dimensional Euclidean space, £ 3 , and using the Gauss divergence theorem, it follows that Eq. (1) is equivalent to the 
integral conservation law: 

£l m -ds = 0 , m = 1,2, 3, 4, (2) 

where S denotes the surface around a space-time control volume (CV) in £3 (e.g., the prism ABC — A ! B ! C ! in Fig. lb) 
and \ m = (F m ,G m ,U m ). To evaluate the surface integrals in Eq. (2), gaussian quadratures are required. For second 
order accuracy, the only gaussian point is the centroid of each space-time hyper-surface (e.g., M in Fig. lb) and the 
integrands are approximated by linear functions. Details are given below: 

The inviscid part of the fluxes F^ +1//2 , G £ +1//2 are generated via a Riemann solver (RS, e.g., Roe’s RS). At time level 
t = nAt , the conservative flow variable vector U" is given at all triangular cell centers (e.g. 0,D,E,F in Fig. la). The 
spatial gradients U x and U y can then be reconstructed via a multidimensional limiter (as described below). The left 
and right ( L and R in Fig. la) states at the center of each surface of the space- time CV are established by a linear Taylor 
expansion from their corresponding cell centers. For instance, at the center, M, of the surface ABB ! A ! (Fig. lb), the R 
and L states are respectively: 

u " +1/2 = V n D + (U«) D Ax + (U>A y + (U")d y , (3) 

U " +1/2 = U£+ (U«) 0 Ax+ (U")oAy+ (U")oy. (4) 

Here, the subscripts O and D indicate which cell center flow variables are used, Ax and Ay correspond to either DM or 
OM. The time derivative, U*, is evaluated at time step n with the Cauchy-Kowalewski/Lax-Wendroff time stepping as 
described in [1]. The spatial gradients, U" and U", are evaluated via a multidimensional limiter as described below in 
§ILC. 

The simple quadrature of the surface integrals in Eq. (2) leads to the update formulation of U at O from time step 

n to n + 1 : 

U ” +1 =U "-X;i K +1 %x)k+Gl +1/2 (ny) k ]Al k , (5) 

where F^, G^, ( n x )k , ( n y )k , k = 1,2,3, are respectively the two flux vectors, and the out-going unit normal vector 
components at the centers of the three cylindrical surfaces ABB f A\ BCC f B f and CAA'C'. The lengths of the edges of 
the triangular cell A ABC, AB , BC and CA, are represented by A /^, k = 1 , 2, 3; As the area of AABC and At the time step 
size. 

In order to suppress the undesired pathological behaviors and stabilize the scheme, a multidimensional dissipation 
term [1, 15] is introduced to the basic scheme (5), forming the enhanced time-accurate upwind (ETAU) scheme. The 
multidimensional dissipation is sketched in §II.D. 

C. Multidimensional limiting for spatial gradients 

In the reconstruction stage, for second order accuracy, U must be reconstructed as a linear vector function within the 
current cell via a limiter. Here, a multidimensional as opposed to dimension by dimension limiter is employed [9, 
10]. Let u denote a component of U, by simple finite difference, three linear equations for the gradients ( u x , u y ) are 
obtained: 

(*£> ~x 0 )u x + (y D - yo)u y = u D - u 0 , 

(xe ~ x 0 )u x + {y E - yo)u y = u E - Mo, 

(xf -*o)M;c + iyF -yo)u y = Up - U 0 . 

Here, the subscripts O and D,E,F denote the corresponding cell centers (Fig. la). Generally, any combination of two 
of these equations yields a set of spatial gradients for m. There are three sets of such gradients, namely, («j^ , uf), i = 

1,2,3. Their moduli or norms, m* = (m ^) 2 + (m}^) 2 , i = 1 , 2, 3, are used as a measure. Then, a multidimensional 

limiter based on the magnitude of the norms is applied to all spatial gradients to achieve a single, unified set of 
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gradients for U. Two such multidimensional limiters are recommended here: the minmod limiter and the extended van 
Albada limiter [16] (weighted averaging). 

For the minmod limiter, the gradient ,u$) corresponding to m* = min(mi,m 2 ,m^) is chosen. When the ex- 
tended van Albada limiter is used, the final gradients are evaluated through weighted averaging, 

(1) , (2) . (3) (1) , (2) , (3) 

_ W\U> X + W 2 U x +TV3»j _ w l u y + w 2 Uy + W 3 U £ 

W\ + W2 + W 3 ’ y W\ + W2 + W 3 

with wi = (w 2 W 3 ) a , W 2 = (wiW 3 ) a ,W 3 = (mi/W 2 ) a , a > 0. The extended van Albada limiter or weighted averaging 
was previously used in the CE/SE method [9, 10]. Different values of a provide numerical dissipation at different 
levels. A small a, e.g.,a = 0.5, usually yields less dissipation, and is appropriate for aeroacoustics computations. In 
the presence of a shock or a contact discontinuity, a larger a(= 2.0) is needed. In the rare extreme situation with high 
Mach number and strong shocks, the computed gradients u x , u y need to be further limited by a factor in order to keep 
the numerical procedure stable. 

After the gradients for each component of U are evaluated, the reconstructed U is linear within the entire grid cell, 
A ABC (Fig. la and b). The solution can be considered as either second order accurate or high resolution. 

D. A multidimensional dissipation model 

Despite their great success, upwind schemes are still prone to suffer from local numerical instability, or the so-called 
‘pathological behavior,’ such as the carbuncle phenomenon, the expansion shock, etc., when the Mach number is high. 
In [1, 15], a multidimensional dissipation model was proposed to systematically suppress such pathological behaviors 
and stabilize the scheme. 

To begin, U = XJ avg can be obtained by averaging U* at time step n at the cell interface centers (mid-points) M, N, P. 
They are the R states for RS extrapolated from their corresponding neighboring cell centers (Fig. lc): 

U=(U^+U^ + U^)/3. 

Then, the multidimensional dissipation is introduced by simply replacing U" with a weighted average ptf + (1 — p)U n 
where p, 0 < p < 1, is the weighing factor for dissipation control. A larger p produces higher numerical dissipation. 
With the dissipation model, Eq. (5) is modified to : 

U” +1 =pu+(l — P)U”— ^j^K +l/2 (n x ) k + G n k +1/2 (n y ) k ]Al k , (6) 

Most subsonic or low supersonic flows are insensitive and resilient to the value of p, and P = 10~ 3 is a reasonable 
choice. For flows with extremely high Mach number and strong shocks, P may need to increase up to 0.3 or higher. 
However, to avoid excessive dissipation, one should keep P as small as possible. 

III. Numerical Examples for Flows at Low Mach Number 

The present ETAU scheme is based on an unstructured grid topology and hence is flexible for complex geometries. 
The scheme covers a broad spectrum of flow from hypersonic to low subsonic, including aeroacoustic computations. 
Two examples with medium and low Mach numbers, of 0.2 and 0.02 are given below. 

A. Flow past a circular cylinder at Mach number of 0.2 

Consider a circular cylinder subject to a flow with a free stream Mach number of 0.2. For convenience, the cylinder 
diameter, the free stream speed of sound, and density are used as scales for nondimensionalization. The Strouhal 
number (nondimensional frequency) is here conveniently defined based on the diameter and the free stream speed of 
sound (Note: usually, a Strouhal number is based on the free stream flow speed). The computational domain contains 
about 91 , 000 triangular cells, spanning between —2.5 < x < 10, and — 5 < y < 5. With van Albada limiter (a = 0.5) 
and P = 0.0001, an initial marching of 80, 000 time steps were performed to ensure the start-up transient has exited the 
domain. In the next 170, 000 time steps the pressure history is recorded at an observation point (— 1 , 2) in the field. The 
point is located outside the vortex streets to ensure that only acoustic wave fluctuations are recorded. Figure 2 shows 
the instantaneous isobars at this final time step. Vortex streets downstream of the cylinder are clearly observed. Then, 
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Figure 2. Mach number 0.2 how past a circular cylinder; left: the unstructured computational domain with unstructured grid , right: snapshot of 
isobars at time step 250,000. 


a fast Fourier transform (FFT) is applied to the pressure history. The results are presented as a PSD (power spectral 
density) plot in Fig. 3. The pressure history consists of 17,000 sampling points (recorded every 10 time steps), which 
is just enough to cover 2 14 = 16, 384 points for the FFT. In Fig. 3, the single peak of the PSD corresponds to a Strouhal 
number, St = 0.037, which compares very well with the experimental data [17] (Table 1). Figure 3 also shows that the 
value of (3 has little influence on the Strouhal numbers at the curve peaks. However, the larger value of p, e.g., 0.3, 
represents larger dissipation. 


PSD at (-1,2), beta=0.0001 


PSD at (-1,2), beta=0.3 



Strwhal No. Strouhal No. 


Figure 3. PSD at a field point (—1,2) with different ps; left: p = 0.0001; right: p = 0.3; showing that p has no influence on the PSD and Strouhal 
No. at the peak (PSD=—436, St.= 0.037). Due to the larger p, the PSD curve on the right decays faster with increasing Strouhal No. 


B. Flow past circular cylinder at a low Mach number of 0.02 

In this case, the computational domain, the grid, and the boundary conditions are identical to the one above in §111. A 
except the Mach number is reduced to 0.02 and the Reynolds number is set at 15, 000. The same parameters, a = 0.5, 
P = 0.001, and At = 0.002 are used. However, as the flow is now much slower, it takes more than 1.5 million time 
steps for the start-up transient flow to pass through the domain. The pressure history at 3 locations is recorded with 
additional 4.1 million time steps. To ensure numerical accuracy for long run times, the computation is carried out in 
double-precision mode. Figure 4 illustrates snapshots of the vortex streets in the wake of the circular cylinder flow 
at time step 2.3 million and 5.7 million. A fast Fourier transform (FFT) is again applied to the pressure history. The 
result is presented as a PSD plot in Fig. 5. There are 16,400 sampling points (recorded every 250 time steps) in the 
pressure history. 
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Figure 4. A Mach number 0.02 flow past a circular cylinder; left : snapshot of vortex streets at time step 2.3 million, right: snapshot of vortex 
streets at time step 5.7 million. 


PSD at (-0.9 ,-2. 5) 



Strouhal No. 


Figure 5. PSD at a held point (—0.9, —2.5), showing a Strouhal No. st = 0.0039, which agrees with the experimental data in [18]. 
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We note that for both test cases, the Strouhal number based on the free stream speed of sound is used, which is 
convenient in numerical computations. In Table 1, for comparison, it is converted back to a standard Strouhal number 
that is based on the free stream flow velocity. Table 1 demonstrates that for Mach numbers of either 0.2 or 0.02, the 
computed frequencies (Strouhal numbers) agree well with the experimental data [17, 18]. 


Table 1. Comparison of experimental and computed aeolian noise frequencies 


Mach No. 

Reynolds No. 

Exp. Freq. 
(Strouhal No.) 

Comput. Freq. 
(Strouhal No.) 

0.2 

90,000 

0.1846 [17] 

0.185 

0.02 

15,000 

0.194 [18] 

0.195 


IV. Examples for Flows at High Mach Number 

The upwind schemes work well with many benchmark type problems at high Mach numbers. However, due to 
some inherent local instability, they may suffer from the general symptom of the so-called pathological behaviors 
described by Quirk [11] and others. The ETAU scheme is designed to systematically suppress such pathological 
behaviors by adding the aforementioned multidimensional numerical dissipation. Various examples of curing the 
pathological behaviors have been demonstrated in [1], rendering the ETAU scheme a simple but robust way for com- 
putation in high Mach number flow regime. The following are two typical examples at high Mach number, namely, 
supersonic flow past circular blunt body and a strong shock diffracting over a sharp edge. 

A. Supersonic fbw past a circular blunt body 

Consider a supersonic flow of Mach number 10 past a circular cylinder (blunt body). A slip condition is imposed 
on the solid wall. Simple extrapolation is applied to the outflow boundaries. A free stream of Mach number 10 is 
imposed at the inflow boundary. For flows with shocks, the parameter, a, in the van Albada limiter is usually set to 2. 
For high Mach number flows with shocks, the dissipation parameter p may be increased from the level of 1 0 ~ 3 to 1 0 -2 
or even 10 -1 . To avoid excessive damping, one should keep P as low as possible. Also for high Mach number flows 
with shocks, we note that for convergence, the spatial gradients of the (conservative) flow variables after applying the 
limiting procedure may need to be further reduced by a factor of c, 0 < c < 1 (see e.g., [19]). 

This problem is tested on both structured and unstructured grids. The first grid presented is actually a triangulated 
structured grid but treated as unstructured data. There are about 80, 000 triangular cells in the computational domain. 
As shown in Fig. 6 a and b, the carbuncle is prominent with p = 0 but disappears when a maximum p =0.05 is applied. 
Here, P is at its maximum value near the stagnation area where the carbuncle phenomenon occurs and gradually tapers 
down to 0.01. 

Recently, Kitamura et al [20] reported that with an upwind scheme, the carbuncle depends on a variety of issues 
such as the free stream Mach number, grid shape, etc. Also, curing of the carbuncle is very difficult on a truly 
unstructured grid. So, our next test is on a truly unstructured grid (Fig. 7). There are about 75, 000 triangular cells in 
the grid. All of the flow conditions remain identical to the case with a triangulated structured grid. 

For the present ETAU scheme, in order to achieve a correct solution as compared to Fig. 6, it is only necessary to 
adjust the maximum value of P . The treatment is simple and effective. When the maximum value of P in the stagnation 
region is increased to 0.25, the carbuncle is completely suppressed. With increased dissipation, the bow shock in the 
stagnation region appears thicker (Fig. 7), which is a trade-off for using the unstructured grid. 

B. Supersonic few past a sharp edge 

Another typical pathological behavior is the expansion shock that forms in supersonic flow of high Mach number past 
a sharp edge. Consider the problem of a strong shock diffracting over a 90° edge [11]. Here, the shock Mach number 
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Figure 6. Supersonic flow of Mach number 10 past a circular blunt body, (a) carbuncle phenomenon with Godunov scheme (P = 0) and a 
triangulated structured grid; (b) with max. p — 0.05 imposed, the carbuncle disappears. 



Figure 7. a: isobars with truly unstructured grid, max. p — 0.25; b: enlargement of the stagnation region with the grid. 
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is 5.09. There are about 14,400 triangular cells in the domain. The flow is initialized to the ambient condition: 

(po,Ko,vo,/?o) = (1,0,0, 1/y), y= 1.4, 
and conditions at a shock Mach number of 5.09 are imposed at the inlet: 

(p uUhVupi) = (5.0294,4.0779,0,21.4710). 

At the top, bottom, and along the surfaces of the rectangular block, a slip wall condition is imposed. A simple 
extrapolation condition is applied to the outlet boundary. Figure 8 shows the density contours with and without the 
cure. With p changing slightly from 0 to 0.001, the expansion shock disappears. Details of the expansion shock and 
expansion fan are shown in Fig. 9. 


a d 




km 

mm 


Figure 8. Supersonic flow diffracted over an edge of 90°; a: (3=0, expansion shock is present, b: p = 0.001, expansion shock disappears. 



Figure 9. Enlargements of the density contours near the edge, showing details of the expansion shock and expansion fan. 


V. Concluding Remarks 

The ETAU scheme [1] was originally designed to suppress the pathological behaviors that haunt the upwind 
schemes. The scheme features a built-in multidimensional dissipation model and adopts concepts from high order 
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upwind schemes, such as accurate evaluation of surface fluxes, multidimensional limiter, etc. The ETAU scheme is 
nominally second order accurate in space and time. Also, employment of the unstructured grids enables flexibility in 
geometry. 

In the present paper, in addition to capability demonstrated in [1], the ETAU scheme also exhibits good behavior 
in a wide range of Mach number, and can be considered an ‘all-speed’ or ‘Mach number insensitive’ scheme. As 
demonstrated in § III and IV, for Mach number from 0.02 to 10 in these examples, the scheme works well, and one 
only needs to adjust one or two parameters, namely, (3 and a. 

Selection of a and P is based on the flow Mach number range, rather than on each individual flow. For flows at 
high Mach numbers, a is usually set to 2.0, otherwise, a should be smaller. Generally, P = 0.001 is good enough for 
most flows. However, for flows with high Mach number and shocks, P may be increased to the level of 10 _1 . To avoid 
excessive dissipation, p should be kept as small as possible. 

With appropriate values of P and a, the pathological behaviors and troubles with truly unstructured grid can be 
removed, while invoking no preconditioning at low Mach numbers. The ETAU scheme is thus rendered as a viable 
(nearly) ‘all-speed’ upwind scheme for both CFD and CAA. 
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